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SUMMARY 

In this paper the effect of thermal escape on the distribution 
of neutral hydrogen in the upper thermosphere is discussed. A 
simplified model is proposed which assumes an extended transition 
region between the collision dominated lower thermosphere and the 
effectively collisionless exosphere. In this region collisions with 
the predominant constituent, oxygen, are still too frequent to be 
neglected, while the effect of thermal escape is strong enough to 
impede the approach to local thermal equilibrium. The behavior 
of hydrogen in the transition region is investigated for different 
oxygen temperatures by means of the Monte Carlo method. Re- 
sults show that stationary conditions are attained in relatively short 
times which vary from about one hour at 1500°K to about one-half 
hour at 2500°K. The velocity distribution of hydrogen is seriously 
perturbed, and shows marked anisotropy as a result of escape; 
consequently, the loss rate of neutral hydrogen from the terres- 
trial atmosphere may be reduced by a factor of three or more in 
the range of temperatures considered. 
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INTRODUCTION 

With the realization that a real gas at finite temperatures is actually an assembly of a vast 
number of particles moving at random velocities, it was recognized that planetary atmospheres 
may be subject to selective thermal escape. This mechanism favors the light elements on account 
of their higher thermal velocities. However, since the early inception of this problem the quanti- 
tative treatment of the effect has not advanced greatly beyond the original crude estimates of 
atmospheric escape rates (References 1, 2, 3, and 4). 

The exploration of space in recent years has renewed and expanded the interest in the subject 
of atmospheric escape. Proper understanding of either the processes involved in the formation of 
the earth’s radiation belts, or a clearer insight into the phenomenon of night- sky Lyman- a radia- 
tion, requires a more precise knowledge of conditions in the earth’s exosphere. Beyond an altitude 
of roughly 1000 km (depending on thermal conditions), the primary constituents of the exosphere 
are now assumed to be the light elements helium and hydrogen. These are precisely the elements 
whose distribution may be subject to effects of escape from the earth’s gravitational field. It is 
the purpose of the present paper to investigate in greater detail than attempted hitherto the effects 
on the distribution of hydrogen just below, and at the base of, the exosphere, and to present results 
which, hopefully, are closer to the actual physical situation than those now widely used in applica- 
tions to problems of atmospheric physics. 

As we see it, the core of the problem lies in the non- equilibrium situation with which escape 
confronts us. As we shall demonstrate in the next section, the well known escape formula of Jeans 
(Reference 4) never holds rigorously, and the effects of the inconsistency appear to have been 
underestimated in the past. At best, escape rates based on this formula represent upper limits to 
the actual loss, while quantities which are sensitive to the actual velocity distribution of hydrogen 
at the base of the exosphere may be more seriously affected. A typical example is the concentration 
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of hydrogen at large geocentric distances (Reference 5). Whereas previous attempts at solving 
the problem of the hydrogen velocity distribution near the critical level have concentrated on the 
high energy tail (References 6, 7, 8, and 9), our treatment covers all velocity ranges (Reference 10). 

For the purpose of this investigation a simple model (suggested by the physical situation) was 
constructed of a semi- infinite container with a semi-permeable upper boundary. Within the frame- 
work of this model, the problem was treated by means of the Monte Carlo method, which appears 
appropriate in view of its basic stochastic character. 

As a first step in attacking the problem, we considered the establishment of a quasi- steady 
state in a layer a few mean free paths in depth extending below the assumed base of the exosphere; 
that is, we assumed the temperature of the oxygen environment, in which hydrogen is embedded at 
such altitudes, to be held constant long enough for hydrogen to approach a steady state and obey 
the equation of continuity. Strong emphasis was placed on consistency of the results with the 
underlying assumptions. 

Results of our computations lend support to the concept of a moderately sharp transition to an 
essentially collisionless exosphere. There is, indeed, a critical region in the atmosphere, just 
below the assumed base of the exosphere, where departure from a collision-dominated situation 
may be observed. It is, however, precisely in this region, where effects of escape compete with 
collisions in shaping the distribution of hydrogen. Furthermore, the fact should not be overlooked 
that the altitude of this critical level shifts considerably as a result of the thermal expansion and 
contraction of the main atmosphere. Consequently, the lower boundary of the exosphere will depart 
to some extent from spherical symmetry. 

Results also bear out the earlier estimate (Reference 10) that the escape rate will be reduced 
considerably, and that the whole velocity distribution will depart from the equilibrium distribution 
appropriate to the ambient atmospheric temperature. The process ensuing from escape may best 
be described as anisotropic cooling of hydrogen, which entails both the anisotropy and reduced 
dispersion of velocities resulting from the rapid loss of fast particles. 

Let us now proceed to a critical review of the conventional approach and search for new 
avenues to a better understanding, suggested by our examination. 

THE CONVENTIONAL APPROACH TO THE PROBLEM AND 
A NEW SIMPLIFIED MODEL 

The escape parameter E at distance R may be expressed by 

mv m 2 
E = 2kT 

where 

R = geocentric distance of the escape level (cm), 
m = mass of hydrogen atom (grams), 
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v ro = escape velocity of a hydrogen atom (cm/sec), 
k = Boltzmann f s constant (1.38 • 10" 23 joules/°K), 

T = temperature at escape level (°K). 

The Jeans escape formula may then be written in terms of the escape parameter, E, as 

F es = <t* 2 N 0 (-^r) 2 (1 + E) e-* (1) 


where 

N 0 = hydrogen concentration at escape level (atom/cm 3 ). 

The Jeans escape formula is obtained by multiplying the velocity distribution, 


f(R, v,0,<£) 



e -mv 2 /2kT v 2 


sin 6 , 


( 2 ) 


by vcostfand integrating over v ro <v <°°, 0 <0 <rr/2, 0<4><2n , and over the spherical surface at 
R as well, where 

v = velocity of hydrogen atom 

and 


0 and 4> - polar angles in velocity space. 


But Equation 1 cannot be exact from a physical point of view because an escape flux implies a 
net flux, and the gas possesses a flow velocity w in the radial direction (Reference 11). Even if 
w were uniform, the velocity distribution in the frame of the planet would obviously be proportional 
to 


exp 


m( v - w ) J 
2kT 


= exp 


_ j 


i(v 2 + w 2 + 2vw cos o) 
2kT 


( 3 ) 


and therefore isotropy is destroyed. While in this frame the above limits of integration are the 
correct ones, the proper results of integration would differ from Equation 1, in view of Equation 3. 

If, on the other hand, we chose to integrate in the moving frame, as is effectively done by putting 
Equation 2 in the integrand, we should have multiplied the velocity distribution f by (v cos 0) -w, 
and changed the limits of integration, both over v and e f from v to (v - w) 2 and from tt/2 to 
cos" 1 (“w/v), respectively. However, Equation 1 may be considered a fair approximation pro- 
vided w « v (implying e » l), and provided v 2 e’^ 2 tends to go to zero rapidly enough. This im- 
plies a highly peaked velocity distribution (i.e., a low temperature for a Maxwellian) which 
contradicts somewhat the requirement that w « v for most values of v. This appears to be one of 
the reasons for the ready acceptance of Jeans' formula, at least by workers in terrestrial aeronomy. 
Until quite recently it was not realized that thermospheric temperatures may go as high as 2000- 
2500 °K, or even higher for short periods of time. 
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Acceptance of Equation 1 implies, moreover, another assumption about the form of the veloc- 
ity distribution at the base of the exosphere, namely that 

f = 0 where < v cos 6 < -v^ 


and 


(4) 


f = fm where -v ro < v cos 6 < 

where f m is the Maxwellian distribution, Equation 2. Then, however, the normalization constant 
used (m/277'kT) 3/ 2 , is no longer strictly valid, again with the result that Equation 1 is no longer 
exact. 

We have gone into the detailed formal discussion of the simple Equation 1 for two main rea- 
sons: First it points out the limitations of this formula's applicability, and secondly the conditions 
expressed by Equation 4 may serve as guidelines for further attempts at a solution. Indeed, if 
one accepts the model of a sharp transition to a collisionless exosphere, one may obtain a set of 
plausible boundary conditions for f at r by relaxing the requirement that f be Maxwellian in the 
interval (-v*, v ro ) . As we shall see below, the results of our model computations agree very well 
with these assumptions. 

We shall not go here into a discussion of further difficulties involved in Equation 1 and at- 
tempts to resolve them, as this has been done by one of the authors elsewhere (Reference 12); 
instead, we proceed directly to a brief exposition of our simplified model and the method applied 
in its solution. 

In consequence of the above discussion the pertinent facts in relation to the escape problem 
appear to be that equilibrium is effectively destroyed by thermal evaporation of the escaping 
element and that near -equilibrium theories, such as are implied by the application of Chapman's 
diffusion equation (References 13 and 14), do not yield the desired information about the actual 
velocity distribution of the escaping element close to the base of the exosphere. 

If one considers, however, some time scales which arise naturally in the context of the escape 
problem, a considerable simplification in the treatment of the problem may be achieved. The two 
significant times we have in mind are effusion time, t e , and diffusion time, t d , to be defined as 
follows. 

Let the base of the exosphere be at altitude z 0 , and the diffusion velocity of the escaping ele- 
ment at altitude z' be w(z'). Sincew(z') = dz'/dt, an average flow velocity (directed upwards in 
the case of hydrogen or helium), the mean time of transport for hydrogen (or helium) from any 
arbitrary altitude z < z Q , to the base of the exosphere is given by 


t d ( z ) 




dz * 
w(z') 


(5) 
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Taking the difference between, say t d ( Zl ) and t d (z 2 ) one obtains, in effect, the time of sojourn of 
a diffusing ensemble between two levels and z 2 . Of course, w(z) has to be obtained from a 
solution of the problem, but we may use as an ad hoc approximation a formula such as Equation 12 
in Reference 14. 

The definition of effusion time, t e , is based on the consideration that since in reality escape 
does not take place from a sharply defined layer, it is reasonable to state that the bulk of escaping 
hydrogen (or helium) will be derived from a layer on the order of an atmospheric scale height, 
below the exosphere, and that in the escape region the mean free path is of about equal dimensions. 
Atoms diffusing from below this layer have a high probability of colliding in it and atoms above it 
are likely to have had their last collision there. Consider now the thermal evaporation of such a 
layer in the absence of sources. This process is governed by the rate equation 


where t?c m ~ 2 is the content per cm 2 of such a layer, given by 77 = nh, where n and h are some mean 
values of the concentration and scale height in this layer. 

On the other hand, &n/dt , is the instantaneous escape flux, given (ad hoc) by Equation 1, whence 


t 

e 



( 1 - E) e E . 


( 7 ) 


It appears that the ratio t d /t e may serve as a reasonable criterion in estimating the departure 
of the distribution from equilibrium. Once a set of conditions is reached where this ratio dwindles 
below unity, escape and diffusion will be of comparable significance in shaping the distribution of 
the escaping element. It suggests that we pay special attention to this limited critical region, and 
furthermore implies a subdivision of the upper atmosphere (above the turbopause) into three main 
regions: 

1. The bulk of the thermosphere, where diffusion theory according to Chapman appears en- 
tirely adequate (Reference 14). 

2. The exosphere, where orbital theory (Reference 11), use of Liouville's equation (Reference 
15), and the collisionless Boltzmann equation (Reference 16) yield equivalent results. 

3. Between these regions an extended transition region of a few mean-free-paths depth. Here, 
the approaches used in the other regions fail. 

From the nature of the criteria chosen it should be evident that the boundaries of the transition 
region are not sharply defined. But in analogy with the theory of the exosphere it appears reason- 
able to postulate sharp boundaries, provided sufficient care is exercised in their choice. Our model 
will be plausible and internally consistent if, within certain limits, the results are not sensitive to 
the particular boundaries chosen. This assures us that the characteristic structure of our results 
is not a spurious boundary effect, but results from intrinsic properties of the physical processes 
investigated in our model. How this can be achieved with confidence in our results will be shown 
later. 
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The fact that all throughout the transition region hydrogen is a minor constituent with a rela- 
tive concentration of less than 3 x 10“ 3 justifies further essential simplification in our model 
(i.e., the neglect of hydrogen-hydrogen interaction). Changes in the hydrogen distribution evolve 
from interaction with the surrounding oxygen medium, which constitutes an essentially undisturbed 
thermal bath. 

A crucial simplification results from a consideration of the trajectories followed by atoms in 
the exosphere: those with sufficient velocity are lost "forever”; those below escape energy, how- 
ever, describe ballistic orbits. In spherical symmetry, we may then assume that for each particle 
in this class leaving the escape level, a corresponding mirror particle enters from the exosphere 
after a certain time delay. This suggests that the effect of ballistic orbits may be simulated by 
the boundary conditions at the top: the escape level may be regarded as a selectively permeable 
membrane . Particles above escape energy pass through it unhindered, while slower atoms may 
be considered to be specularly reflected with an appropriate time delay. Once a steady state has 
been approached, of course, there is no time delay and outgoing atoms are exactly balanced by 
incoming particles. 

This property becomes especially simple, mathematically, when in view of the small thickness 
of the transition region a plane-parallel approximation is adopted. The nature of the boundary 
conditions at the top is very similar to those of Equation 4. The main distinction, however, is that 
no particular form of the distribution function can be postulated, as was done there. All that one 
can observe is that near the boundary the distribution function should be symmetrical in the ver- 
tical component of velocity, v 2 , in the domain (-v OT , v ro ) . 

Finally the fact that particles are fed into the region of interest from below by diffusion can 
be taken into account by postulating a source, located at the lower boundary, and (in the steady 
state case) injecting particles at a constant rate. The actual form of the source flux is of no great 
import, however, since in the lower portion of the transition region collisions are still frequent 
enough to bring about a fair amount of randomization. Our results bear out the contention that 
apart from a region in the vicinity of the source, the steady state distribution is not sensitive to 
the particular form of the source flux. 

The considerations expressed suggest the simple physical model we were searching for: a 
container, extending infinitely in the x- y directions and of finite height in z. The container is filled 
with a major component gas in thermal and diffusive equilibrium, subject to a spatial distribution 
determined by gravity and temperature alone. After injection by a source at the bottom, atoms of 
the minor constituent percolate through this medium until they reach the semipermeable top of the 
container where (depending on their velocity) they may be either reflected or pass through. 

SOLUTION OF THE TRANSPORT PROBLEM BY MEANS OF 
THE MONTE CARLO METHOD 

Though the model just delineated appears rather simple and well defined, solution of the cor- 
responding transport equation (the linear Boltzmann equation) by analytical methods is most 
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difficult without surrendering much of its physical content by simplifying approximations. The 
main difficulty lies in the highly singular boundary conditions such as the tripartite division of the 
velocity range in the vicinity of the escape level (-oo ? - v oo) ? (-vco, v<»), (v<», °°), and the finite extent 
of the vertical dimension. The small depth of the container (representing the transition region) 
indicates, however, that the Monte Carlo method may be applied with a good chance of success. 

A previous Monte Carlo treatment of hydrogen diffusion (Reference 17) did not take into ac- 
count escape, and, with its implicit use of the Maxwell- Boltzmann distribution for hydrogen, did 
not aim at elucidating the question of the velocity distribution near the escape level. In distinction 
to that treatment which involved an optical model of the atmosphere with varying transmissivity 
for hydrogen, we apply the Monte Carlo method in a straightforward fashion: we follow a large 
number of particles through their life history in the transition region. The particles are repre- 
sented by their parameters of interest such as speed, direction, position and so on. While the 
particles proceed through the region, relevant parameters are recorded intermittently. When a 
sufficient number of particles have been processed, an appropriate census is taken yielding what 
are essentially histograms representing the distribution of the sample particles in space, velocity, 
direction and so on. A good general account of the technique was given by Cashwell and Everett 
(Reference 18) which obviates a detailed account here. 

Here we wish only to outline briefly the procedure followed in our Monte Carlo calculation. 

In accordance with the assumptions made previously, the hydrogen problem is linear (i.e., each 
hydrogen atom moves independently) interacting only with the particles of the medium. This per- 
mits us to follow each atones life history independently of the others. 

We inject hydrogen atoms at the source (z = 0) at a random time between 0 and a fixed time r m 
The source parameters are drawn randomly from an assumed source parameter distribution. The 
source parameters are scalar velocity, zenith angle, and ’’statistical weight”. The meaning and 
significance of the latter will be explained further on. Next, dependent upon our assumptions with 
respect to the kinetic cross section of hydrogen, we draw a random free path to determine the 
position of the next encounter with an oxygen atom from a free path distribution appropriate for 
the medium. If the position so determined is outside the "container” the life history of the particle 
is terminated without investigating the effects of the projected encounter. If, however, the collision 
is valid, the parameters of the oxygen atom serving as the hydrogen’s partner in the collision are 
drawn at random consistent with the assumed oxygen equilibrium distribution. The change in the 
hydrogen atom’s parameters is determined from the mechanics of the collision (here, we assume 
isotropic scattering of the hydrogen in the center of mass system). A new free path is drawn and 
the procedure is repeated until the H-atom finally leaves the system by either escape at the top or 
diffusion through the bottom of the container. 

The container is divided into layers of width A z, each labeled z if the velocity range is divided 
into intervals v. , and intervals are established for the range of zenith angles also. At fixed times 
(multiples of the fundamental time unit, r) the computation is interrupted and the parameters of 
the H-atom are recorded. A number of "counter arrays" is set up (each counter is empty at the 
beginning of the computations, and corresponds to some combination of intervals of space, velocity 
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and time previously mentioned), and at recording time t k the statistical weight of the particle is 
added to the contents of the counter which correspond to the particle T s parameters (e.g., if n(v i , 
z., t k ) denotes the weight of particles with velocity in v. in the layer z. at time t k , a particle's 
weight will be added to it, provided that the particle at time, t k , possesses the corresponding 
variables of motion). Because of the assumed linearity of the process, the arrays play a passive 
role and do not affect the subsequent motion of particles. 

The stationarity of the process (i.e., our assumption about the invariability of the oxygen 
medium) permits a further essential simplification since it enables us not only to superpose the 
records of different atoms (from linearity), but to superpose records for different times also. 

This is best explained in the language of stochastic processes. 

Let n i (k)be the number of particles with parameters of motion, symbolized by i (at time 
t = t k ) be divided by the total number of particles in the system. Then, in the language of Markov 
processes, it may be termed an "absolute probability" (Reference 19) given by 


n i ( k ) - 2^ n i (°)Pj C i k) (8) 

j 

where n j (0)is the initial distribution, and Pj<. k > are the stationary, k-step, transition probabilities 
from state j to state i . 

Consider now the sum 


Pi (m) - (m+ir'^n^k) = (0) Pj C i k J 


(9) 


k = 0 


i k 


It is easily shown that 0 <P i <1, and that Z] P. = l, so that p. satisfies the conditions for a 
probability. 

Now let N. (m)be the number of particles in state j at timem. In the framework of our model, 
this is derived from contributions of the source at times 0, 1, • • ■ ,m - 1. It may then be repres- 
ented in the form 


n j < m ) ~ ^ n i (°) (pi ( j m) ± * ■ * + p/j 2) + PiV } ) 


tn 

■ ll n i (°> Pi f j n> • 


( 10 ) 
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Nj (m) is, therefore, of the form of Pj (apart from a constant), and instead of adding the contribu- 
tions from the invariable source , characterized by the distribution rc. (0 ) , we may sum over the 
number of particles in state j, at various times up to m, as is done in Equation 9. 

Note that N^iois a monotonous, non-decreasing function (disregarding fluctuations), but these 
are of no interest to us at the moment. Now, if we had an infinite number of states and an infinite 
number of particles, the distribution at m would yield the required stationary distribution. In prac- 
tice, we process a finite number of particles, and the lifetime of particles in the container will 
have a finite maximum, T m . The actual value of T m , obtained in the Monte Carlo computation will 
depend on the distribution of lifetimes and, therefore, will be a function of the number of particles 
processed, as will be the values for the various p.< n >, which we implicitly apply in the process. 
Various criteria may serve us to determine the sufficiency of the sample processed, and we shall 
return to this problem farther on. If we assume that a sufficiently large sample has been processed, 
we observe from Equation 10 that an approximately stationary state has been reached. After T m the 
distribution will not change since no further contributions from the initial distribution are obtained. 

A further excursion into probability theory will clarify our use of the statistical weight and 
the assumed source distribution. Under stationary conditions, we might consider escape at the 
top, not as a final exit from the system but implying transition from there to the source because 
for each particle escaping another particle has to diffuse into the system. When seen in this light, 
and in view of the large but finite number of states considered, our model simulates approximately 
a finite, irreducible, aperiodic Markov chain, which by a theorem of Markov (Reference 20) is 
ergodic. To test the independence of the stationary distribution of the initial distribution (i.e., the 
source), we use effectively two different source distributions by applying the following ruse: we 
draw each particle’s parameters from a source flux distribution corresponding to a Maxwellian at 
the equilibrium temperature, but assign to the particle both a statistical weight 1, as well as a dif- 
ferent statistical weight W which, in effect, alters the source flux distribution in any desired way. 

As we shall show, the distribution of hydrogen in the transition region is practically independent 
of the source distribution. 

A few words may also be appropriate at this point about both the free path routine used in our 
computations and the kinetic cross sections utilized in the routine. To achieve a more realistic 
representation we have taken into account the dependence of cross section on velocity by applying 
some results obtained in the investigation of transport properties of real 1 gases in recent years by 
Mason and co-workers (Reference 21). Among others, collision integrals arising in this context 
have been computed by Krupenie et al (Reference 22) for O-H interactions, and these may be 
readily converted into velocity- dependent cross section (Reference 12). This reduction leads to 
the approximate form <j = 8. lx 1(T 16 x v ~ 0 * 6 cm 2 , where v is measured in kilometers per second. 

It should be noted that by assuming isotropic scattering in the center of mass frame we have 
actually represented an atom as a hard sphere, whose diameter is a function of velocity. 

We assume that the oxygen in the transition region is in equilibrium at constant temperature. 
Throughout this region, the extent of which is small in comparison with the earth’s radius, gravi- 
tational acceleration varies only slightly. One is then justified to assume an oxygen concentration 
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varying exponentially with constant scale height, H, such that, 


N o(z) - -N 0 e -jf (11) 

where N 0 denotes the oxygen concentration at the base of the container, located at z = 0, k(z) is 
the mean free path at z for some particular velocity, is the magnitude of the mean free path at 
z = 0, 6 is the zenith angle of a moving particle, and x = cos <9 . If q denotes the survival probabil- 
ity of a particle, that is, the probability of not making a collision over a travel distance the 
change in q over a pathlength, df, is 


d q _ _ T" d£ - _ i^" dz • 

Now, 


( 12 ) 


(X)' 1 = Na = N ae‘‘ /H = \~ 1 e~ z/H 

N ' O U 


(13) 


Then 


which readily integrates to 


dq 

q 


e~ z/H 
K o x 


dz 


(14) 


q(z) - e ”(H/Aox) (l e Z//H ) . 

Obviously, the probability of a collision up to z is given by 


(15) 


P(z) - 1 - q(z) = 1 - e "(H/A 0 x)(l e z/H ) 


(16) 


P(Z) is a random variable with uniform distribution in (0,1), and the principle of Monte Carlo is 
readily applied, yielding for the altitude of the first collision of a particle starting at z = o 


1 r i 

- H In 1 - Ins , 


(17) 


where s is a random number, between 0 and 1, with the proviso that 


s > e 


-hA„ 


(18) 


If the particle is located at, say, z = z', the location of the next collision is given by a formula 
analogous to Equation 17 with K Q replaced by k Q e z/H , so that 


z 


\ 0 X 

H ln 1 - -jj— e l/H In s 


(19) 
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It should be mentioned that in our model calculation we have not taken into account the effect of 
gravity on the curvature of particle trajectories. This does not mean that gravity has been neg- 
lected altogether. The main effect of gravity is the varying concentration of oxygen. In terms of 
the Boltzmann equation one might say that gravity appears in the collision operator, while the 
force term v • df/dv has been deleted. 


RESULTS OF THE MONTE CARLO COMPUTATION 


The Monte Carlo treatment outlined, in contrast to analytical treatments, requires the in- 
troduction of certain physical parameters in the course of the calculation itself. Before presenting 
the results of the computation, we shall describe briefly the atmospheric model used in obtaining 
our results. The crucial parameters are the temperature of the medium (predominantly oxygen) and 
the number density of ambient oxygen. The large amount of information gained in recent years 
permits a high degree of confidence in the agreement of numerical results with the actual phys- 
ical conditions prevailing in the upper atmosphere. 


For computation in this study, two typical cases were chosen: an isothermal atmosphere at 
1500°K, and one at 2500TC. This choice was motivated by the desire to gain insight into the escape 
problem both under average conditions, and under rather extreme conditions at solar maximum 
represented by the higher temperature. The atmospheric data for the 1500° case are based on 
a variety of sources (References 14, 23, and 24) which show satisfactory numerical agreement 
between the various authors. 


Data for the high temperature case 
were deduced by extrapolation from the 
cited study of Harris and Priester (Ref- 
erence 24) of the time-dependent at- 
mosphere. The relevant parameters, 

N 0 , X 0 , and H are represented in Table 1. 

As mentioned previously, the choice 
of proper boundaries for the computa- 
tion is intimately connected with the con- 
sistency of the model. The effect of 
escape on the distribution should not be 
critically dependent on the particular 
choice of its altitude because, as we know, the sharp escape level is a convenient simplification 
and should not be taken too literally. In practice, this will be a rather extended diffuse region. 

If we are able to demonstrate that departure from approximate diffusive equilibrium (ascribed 
to the effect of escape) occurs independently of the exact boundary, our confidence in the meaning- 
fulness of our results will be greatly strengthened. Figures 1 and 2 show the variation of the 
hydrogen concentration for both temperatures and various upper boundaries. The statistics of 
these computations are not too accurate, but apart from curve 1 in Figure 1, they show good 


Table 1 


Atmospheric Parameters Used in the Computation 


Atmospheric 

Temperature 

(K°) 

Altitude 

(km) 

Oxygen 
Concentration 
N 0 (cm' 3 ) 

Mean 
Free Path 
\ 0 (km) 

Scale 
Height 
H (km) 

1500 

350 

4.12k 10 s 

10.11 

77.07 


700 

4 .35 x 10 6 

1049.34 

2500 

575 

2.00 x 10 s 

25.00 

123.30 


1025 

4.94 x 10 6 

1012.00 
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RELATIVE CONCENTRATION 


agreement in the upper part of the transition region. They also display a significant change in 
slope of the concentration curve to occur in a definite region, the location of which is independent 
of the boundaries chosen. This leads to the conclusion that the increasing dilution of the ambient 
oxygen results in a critical region where escape affects significantly the distribution of hydrogen. 
Application of the usual criteria for the choice of the base of the exosphere (References 11 and 12) 
reveals that this region is located close to, but below this level. This is consistent with our ob- 
servation that escape will result in a changed distribution at the exospheric base, which has to be 
taken into account before applying methods appropriate to the collisionless domain. In the low 
temperature case, this break is quite sharp with the critical region located in the vicinity of 500 
km, while the exospheric base, in agreement with our atmospheric model, it is around 580 km. 

The transition in the 2500° case is less pronounced, occurring around 700 km; its relative smooth- 
ness appears to be a consequence of the considerably increased scale height of oxygen. As Fig- 
ure 2 reveals, the choice of the lower boundary has some import on the results in the lower part 
of the region of interest (but does not affect the upper portion). If we wish to join our results to 
a curve obtained from approximate diffusive equilibrium, this choice has to be done with sufficient 
care. We cannot go into the detailed considerations of this problem here, but the choice of the 
particular level adopted has little effect on the velocity distribution since the concentration of the 

predominant constituent in the lower region is 



Figure 1 —Relative hydrogen concentration versus altitude Figure 2— Relative hydrogen concentration versus altitude 

for different boundaries. for different boundaries. 
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In accordance with the above considerations, the levels chosen as boundaries in our computa- 
tion were 350 km and 700 km for the 1500° case, and 575 km and 1025 km for the 2500 case. 

Of great significance, of course, is the problem of sample size required for processing in a 
computation of this character and the statistical accuracy achieved. As we observed earlier, our 
approximation to a realistic picture of a stationary state will depend quite strongly on the size of 
the sample. Continuous, endless replication of the Monte Carlo procedure should yield exact 
results within the framework of the underlying physical assumptions, but the approach to limiting 
accuracy becomes gradually slower and slower. It is necessary, therefore, to choose some suit- 
able criteria for termination of the computations. Since the problem at hand is of greater com- 
plexity than relatively simple cases amenable to well defined mthods such as sequential analysis 
(Reference 25), we were forced to be guided by some intuitively plausible criteria although we 
were unable to put them into rigorous mathematical terms. It has to be borne in mind, moreover, 
that for given sample size, the accuracy achieved will vary with the nature of the statistic inves- 
tigated. The sample size required for a realistic estimate of the relative concentration, for in- 
stance, will be much smaller than that necessary for a good estimate of the angular distribution 
of particles in a particular velocity range and a particular layer. Since we were interested in a 
great range of statistics, our sample size is on the order of hundreds of thousands of source par- 
ticles, the bulk of which does not penetrate too far into the medium, and is swiftly scattered back 
through its lower boundary. 

An idea of the stability and accuracy obtained in the computations can be gained from inspec- 
tion of Tables 2 and 3, which represent the relative concentration in selected layers of the medium 
versus the number of source particles processed. It appears that for this particular statistic we 
are within an accuracy of 1%. 


Table 2 

Variation of the Relative Hydrogen Concentration in Selected Layers 
with the Number of Source Particles in the 1500° Case. 


Number of Source 


Layer (kilometers) 


Particles 

378-385 

448-455 

518-525 

588-595 

658-665 

20036 

.939347 

1.035000 

.998599 

.800252 

.700939 

46173 

.979458 

.990204 

.996831 

.796101 

.694602 

64173 

.969872 

.984525 

.994282 

.802329 

.695022 

90000 

.970037 

.989772 

.999978 

.807954 

.697872 

120000 

.976554 

1.000267 

1.011560 

.814299 

.706973 

150000 

.994485 

1.015202 

1.004433 

.837367 

.723087 

180000 

1.000058 

1.014311 

.997391 

.843658 

.726010 

210000 

1.012421 

1.012463 

.994761 

.843256 

.724136 

240000 

1.010880 

1.010340 1 

.986317 

.846516 

.717696 

270000 

1.006388 

1.000687 | 

.968360 

.834693 

.703233 

300000 

1.015932 

1.008722 | 

.976175 1 

.844235 

.715809 

330000 

1.010353 

1.004137 

.970938 

.841901 

.714477 

360000 

1.012284 

1.007843 

.972897 

.847781 

.719022 

390000 

1.016374 

1.010825 : 

.970580 

.848068 

.716457 
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Table 3 


Variation of the Relative Hydrogen Concentration in Selected Layers 
with the Number of Source Particles in the 2500° Case. 




Layer 

(kilometers) 


Number of Source 
Particles 

611-620 

701-710 

791-800 

881-890 

971-980 

21500 

1.043191 

.917147 

.860685 

.744332 

.650085 

67000 

.963635 

.897185 

.842384 

.747455 

.616731 

100000 

.973088 

.904620 

.858624 

.751596 

.629521 

130000 

.979538 

.913705 

.855260 

.754563 

.632597 

160000 

.979349 

.915424 

.857541 

.754975 

.632721 

190000 

.976984 

.913239 

.856988 

.753267 

.634049 

220000 

.975521 

.910686 

.853568 

.751488 

.633437 

250000 

.974909 

.911627 

.847652 

.746549 

.629097 

280000 

.975148 

.907655 

.845661 

.743627 

.625194 

310000 

.979932 

.911077 

.847217 

.743735 

.626639 

350000 

.973840 

.906308 

.842933 

.741363 

.623740 


The results of the computation are best classified roughly into two categories: microscopic 
results and macroscopic results. By the former we mean, primarily, the detailed velocity distri- 
bution as a function of altitude, i.e., quantities which are not likely to be amenable to direct meas- 
urement in the near future. 

We present first, therefore, the salient features of the macroscopic results such as the varia- 
tion of hydrogen concentration with altitude and typical relaxation times as well as diffusion times 
and velocities. 



Figure 3 — Relative hydrogen concentration versus altitude 
fori 500 °K. 


Figures 3 and 4 present the relative con- 
centration of hydrogen versus altitude in the 
transition region for the two cases treated. 
As can be seen, the departure from diffusive 
equilibrium is not too pronounced. The sig- 
nificance of the departure lies, however, in 
the increased rate of change as one ascends 
in the transition layer. Previous work (Ref- 
erences 14 and 26) implicitly assumed that 
with increasing altitude and growing distance 
of the physical source of hydrogen at about 
100 km, the approach to diffusive equilibrium 
becomes closer and closer. Our results in- 
dicate that this is not the case; equilibrium 
cannot be attained when collisions become in- 
creasingly rare, and the effects of selective 
removal make themselves felt. As Figure 4 
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(t)/nf 


shows the smooth fit to an equilibrium curve 
at the lower end has been achieved here with a 
high degree of success. In the other case we 
were not quite aware of the implications of this 
problem when calculations were started, and in 
this case our results should be fitted to the dif- 
fusion curve at a level between 450 and 500 km. 

As remarked above, this does not significantly 
affect the results from the critical region around 
500 km on. 

Figure 5 presents a typical curve of how 
the concentration in a layer varies with time. 

It is seen that it is very close to a time depend- 
ence of the form 1 - e~ t/r . If, on this basis, we 
compute the numerical values of the relaxation 
times, we obtain from a least square calculation for selected layers the numerical values presented 
in Tables 4 and 5. As is to be expected, the lower layers, where collisions are frequent, approach 
stationary conditions faster than the upper, rarer layers. If we denote by n(t) the number content 
of the transition region at time t , and by n f the content upon reaching a steady state, a plot of 
l - [n(t)/n f ] is of the form presented in Figure 6. As we observe from these curves, the times for 
attainment of a high degree of stationarity, say 90%, vary from~ 2500 seconds at 2500° to^ 1500 
seconds at 2500°. The short times involved indicate that the distribution of the minor constituent 
follows very closely the thermal conditions prevailing in the main atmosphere, which should be 
taken into account when attempting to calculate the lateral flow of hydrogen due to concentration 
gradients (References 27 and 28). In reality, due to the time variations of the main atmosphere 
(which seem to be well understood (References 24 and 29)) stationary conditions will never be 



Figure 4 —Re I ati ve hydrogen concentration versus a I ti tude 
for 2500 °K. 




Figure 5 — Typical curve showing the hydrogen 
concentration-versus-time approach to steady state. 


Figure 6— Transition region-versus-time approach 
to steady state. 
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Table 4 


Relaxation Times of Selected Layers in the 1500° Case. 


Layer (km) 

378-385 

448-455 

518-525 

588-595 

658-665 

t (sec) 

943.0 

1096.0 

1150.0 

1222.4 

1192.1 


Table 5 

Relaxation Times of Selected Layers in the 2500° Case. 


Layer (km) 

611-620 

701-710 

791-800 

881-890 

971-980 

t (sec) 

570.6 

637.6 

660.1 

666.2 

659.2 


fully attained, but the approach to them appears faster than was thought hitherto, which would 
prevent the smoothing out of local differences due to lateral flow. A final answer to this problem 
can be obtained only after extension of this study to the time dependent case. 

In regard to the problem of diffusion times, in our computation we have obtained the mean 
time of travel of escaping particles between the lower and upper boundaries of the transition re- 
gion. These correspond to mean diffusion velocities as follows: 

~3.4 x 10 4 cm sec -1 for 1500° case 

~9.0 x 10 4 cm sec -1 for 2500° case. 

Note should be taken that these mean diffusion velocities are not identical with the mean flow 
velocity of hydrogen, except, ideally at the escape level itself where escaping particles provide 
all the net flow. The former represents an average taken only over those particles which ultimately 
pass through the upper boundary, and does not take into account those hydrogen atoms which are 
turned back at one or another location in the "container" thus terminating their life history by 
passing downward through the lower boundary. 

We now turn to a review of the main microscopic results. In our computations we obtained the 
distributions of various parameters in selected layers spaced about equally throughout the region 
of interest. The width of these layers was chosen such that their kinetic depth varied from about a 
mean free path near the lower boundary to a small fraction of a mean free path near the upper 
boundary; this is in accordance with the assumption that changes in the distribution will be faster 
near the escape level. The parameters whose distribution was obtained, were the following: ver- 
tical velocity, horizontal velocity, total velocity, and angular distribution of various velocity groups. 
For the distribution of the vertical velocity component, the double bookkeeping procedure outlined 
previously was adopted. 

Samples of the distributions are presented in Figures 7 through 12. 
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VERTICAL VELOCITY DISTRIBUTION AT 580 km m VERTICAL VELOCITY DISTRIBUTION AT 350 km 







VELOCITY DISTRIBUTION AT 580 km VELOCITY DISTRIBUTION AT 350 km 



SCALAR VELOCITY ( 1 unit = 4.975 km) 

Figure 1 1 —Distribution of scalar velocity for 1500°K case (dashed curve - Maxwellian at 1500°K). 
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These figures are arranged in pairs; one member of the pair referring to the distribution at 
a level close to the lower boundary, the other referring to a level in the vicinity of the assumed 
level of the exospheric base. In all these figures a comparison curve is included, showing a dis- 
tribution corresponding to equilibrium at the temperature of the ambient oxygen, and normalized 
to the area under the curve obtained in our computations. The unit of velocity used in these curves 
is (2kT/m) 1/2 (i.e., 4.975 km/sec at 1500°K and 6.422 km/sec at 2500°K). Vertical units are arbi- 
trary, but are the same for both members of a .pair so that the ratio of their areas corresponds 
to the ratio of hydrogen concentrations at the corresponding levels. 

Figure 7 shows the computed distribution of the vertical velocity component, v z , around 350 km 
and 590 km for the 1500° case. At the lower level the fit to the assumed Maxwellian is almost 
perfect. In the computed curve, there is only a slight asymmetry which is needed to sustain a net 
flow. At the higher level where both the departure from equilibrium and the asymmetry become 
quite pronounced, the disagreement between the curves mount with increasing velocity. Even more 
striking is this effect in the high temperature case presented in Figure 8, where the depletion of 
the high velocity tail of the distribution can be better observed. 

Figures 9 and 10 display the behavior of the horizontal velocity for both cases treated. Evi- 
dently, in view of the dependence of escape on direction, the distribution of this component is far 
less altitude dependent than that of the vertical component. Corresponding figures for the distribu- 
tion of total velocity are Figures 11 and 12, which reveal a degree of departure from equilibrium 
intermediate between those of the separate distributions. 

In the absence of an equilibrium distribution, isotropy is destroyed and the angular distribution 
becomes a function of velocity. Insight into this effect can be gained from Figure 13 which presents 
the angular distribution of narrow velocity groups centered around v = 1. 1 and v = 2. l (in the units 
described above) for the high temperature case. The lower velocity is near the mean velocity, 
while the higher velocity is beyond escape velocity. The marked change in the angular distribution 
is characteristic of the process studied: in the lower region, for low and moderate velocities, the 
angular distribution is isotropic for all practical purposes; in the higher velocity range a slight 
degree of anisotropy may be observed. In the escape region, anisotropy is destroyed for all veloc- 
ities, yet, in the velocity range below escape velocity, symmetry in the forward and backward direc- 
tions is conserved while beyond this limit asymmetry is most pronounced; this is in agreement 
with the assumptions made earlier. 

To obtain a more quantitative estimate of the departure from equilibrium the results were sub- 
jected to thorough statistical analysis, which confirmed the trend in behavior of the escaping gas. 

An example of this is presented in Table 6. In this table we present the results of a y 2 test, as 
well as a Kolmogorov-Smirnov test for the hypothesis that our computed distribution is a sample 
of an equilibrium distribution at the temperature of the oxygen medium. The numbers appearing 
in the table are, in effect, the probabilities that this hypothesis is true in selected layers of the 
transition region. 
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Figure 13— Typical Angular Distribution Curves. 


Table 6 


Significance Levels for the Distribution of Vertical and Horizontal Velocity at 1500° K. 



Vertical Velocity 

Horizontal Velocity 

Layer (km) 

Y 2 

K-S 

Y 2 

K-S 

350-357 

0.720 

0.500 

0.090 

0.085 

378-385 

0.310 

0.027 

0.035 

0.029 

413-420 

0.360 

0.250 

0.042 

0.051 

448-455 

0.038 

0.025 

0.013 

0.021 

483-490 

0.002 

0.001 

0.032 

0.027 

518-525 

0.001 

<0.001 

0.024 

0.030 

553-560 

<0.001 

<0.001 

0.009 

0.007 

588-595 

<0.001 

<0.001 

0.011 

0.008 

623-630 

« 0.001 

«0.001 

0.010 

0.010 

658-655 

« 0.001 

«0.001 

0.007 

0.005 

693-700 

« 0.001 

«&0.001 

0.006 

0.003 
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Of great value in our interpretation of the escape process, are further results of the statistical 
analysis. Since temperature is, in essence, a measure of the dispersion of velocities, we may go 
one step further and define an effective temperature for the separate directions of velocity as 
measured in terms of the dispersion of these velocity components. The results of this approach 
are presented in Figures 14 and 15. It may be observed from these figures that the vertical ef- 



350 400 500 600 700 


ALTITUDE (km) 

Figure 14— Effective temperatures versus altitude for 
1500°K case. 



Figure 1 5— Effective temperatures versus altitude for 
2500 °K case. 


fective temperature decreases steadily, whereas, 
the horizontal effective temperature decreases 
initially and remains more or less constant 
thereafter. This is what we term anisotropic 
cooling of the escaping gas; in other words, we 
might say that temperature is a function of 
direction and that the rate of change of temper- 
ature varies with direction also. In any further 
investigation of the escape problem this aspect 
may well serve as a convenient starting point. 

As a final sample of our results, we pre- 
sent Figure 16, a graph of the normalized net 
flux in the transition region for the case of 
2500°. As required by steady state conditions, 
this is constant within the limits of our statis- 
tical accuracy; it is another indication of the 
consistency of our results. No such curve is 
presented for the 1500° case, because in view 
of its small magnitude, the statistical accuracy 
in our determination of the mean flow velocity 
is much poorer. 

We close with a remark pertaining to the. 
loss rate of hydrogen. Table 7 presents the 
effusion velocities of hydrogen at the base of 



Figure 16— Normalized net hydrogen flux versus altitude 
for 2500 °K case. 
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the exosphere based on the calculations of 
Kockartsand Nicolet (Reference 14), as well 
as the results of the present calculations. 


Table 7 

Effusion Velocity of Neutral Hydrogen 


Temperature of 

Effusion Velocity (cm/ sec) 


the Ambient 
Atmosphere (°K) 

Kockart & 
Nicolet 

Present 

Calculations 

Ratio 

1500 

7.5x 10 3 

~ 1.0 X 10 3 

7.5 

2500 

4.0 x 10 4 

~ 1.5 X 10 4 

2.67 


It may be somewhat surprising that the 
(relative) discrepancy is greater in the low 
temperature case, but one should bear in 
mind that in this case the quantities involved 
are so small that a small absolute change 
may produce a great relative change. More- 
over, the escape process is, to a certain 
extent, self-governing. The changes in the 

distribution resulting from escape reduce escape. Within the limits of our model's validity, we 
observe that the loss rates based on equilibrium appear to be an overestimate by a factor of 3 
or more in the range of temperatures covered by this study. 
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